#include "stdafx.h"

#include "hnum_pzsp_defs.h"

namespace harlinn
{
    namespace numerics
    {
        namespace SuperLU
        {
            namespace DoubleComplex
            {

                void
                pzgstrf_panel_dfs(
		                  const int  pnum,  /* process number */
		                  const int  m,     /* number of rows in the matrix */
		                  const int  w,     /* current panel width */
		                  const int  jcol,  /* leading column of the current panel */
		                  SuperMatrix *A,   /* original matrix */
		                  int *perm_r, /* row pivotings that are done so far */
		                  int *xprune, /* in */
		                  int *ispruned,   /* in */
		                  int *lbusy,      /* in; size n */
		                  int *nseg,	   /* out */
		                  int *panel_lsub, /* out */
		                  int *w_lsub_end, /* out; values irrelevant on entry */
		                  int *segrep,     /* out */
		                  int *repfnz,     /* out */
		                  int *marker,     /* modified */
		                  int *spa_marker, /* modified; size n-by-w */
		                  int        *parent,     /* working array */
		                  int *xplore,     /* working array */
		                  doublecomplex *dense,      /* out; size n-by-w */
		                  GlobalLU_t *Glu         /* modified */
		                  )
                {
                /*
                 * -- SuperLU MT routine (version 2.0) --
                 * Lawrence Berkeley National Lab, Univ. of California Berkeley,
                 * and Xerox Palo Alto Research Center.
                 * September 10, 2007
                 *
                 * Purpose
                 * =======
                 *
                 *   Performs a symbolic factorization on a panel of columns [jcol, jcol+w).
                 *   It skips all those busy descendants that are worked on by other
                 *   processors along the e-tree path.
                 *
                 * Notes
                 * =====
                 *
                 * (1) panel_lsub[0:w*n-1]: temporary for the nonzero row indices below 
                 *     the panel diagonal, which will be used later in the inner LU
                 *     factorization. For the busy columns, some of the nonzeros in U
                 *     may be mistakenly placed in this list, because "perm_r" is
                 *     still "empty". Later, during dcolumn_dfs in the inner factorization,
                 *     we must filter those nonzeros belonging in U.
                 *
                 * (2) A supernode representative is the last column of a supernode.
                 *     The nonzeros in U[*,j] are segments that end at supernodal
                 *     representatives.
                 *
                 * (3) The routine returns one list of the supernodal representatives
                 *     in topological order of the DFS that generates them. This list is
                 *     a superset of the topological order of each individual column within
                 *     the panel. The location of the first nonzero in each supernodal
                 *     segment (supernodal entry location) is also returned. Each column
                 *     has a separate list for this purpose.
                 *
                 * (4) Two marker arrays are used to facilitate dfs:
                 *     marker[i] == jj, if i was visited during dfs of current column jj;
                 *     marker1[i] == jcol, if i was visited by earlier columns in this panel;
                 *
                 * (5) The dfs stack is the combination of xplore[2*m] and parent[m]:
                 *     xplore[k]     - pointer to k's adjancency list where search begins
                 *     xplore[m + k] - pointer to k's adjancency list where search ends
                 *
                 * (6) Array mappings
                 *     marker: A-row --> A-row/col (0/1)
                 *     repfnz: SuperA-col --> PA-row
                 *     parent: SuperA-col --> SuperA-col
                 *     xplore: SuperA-col --> index to L-structure
                 *
                 */
                    NCPformat *Astore;
                    doublecomplex    *a;
                    int       *asub;
                    int       *xa_begin, *xa_end;
                    register int krep, chperm, chmark, chrep, kchild, myfnz;
                    register int k, krow, kmark, kperm, fsupc;
                    register int xdfs, maxdfs, kpar, jj, nextp;
                    register int nextl_col;/* next open position in panel_lsub[*,jj] */
                    int       *marker1;	   /* marker1[jj] == jcol if vertex jj was visited 
			                      by a previous column within this panel.   */
                    int       *repfnz_col; /* start of each column in the panel */
                    doublecomplex    *dense_col;  /* start of each column in the panel */
                    int       *xsup, *xsup_end, *supno, *lsub, *xlsub, *xlsub_end;

                    int       *col_marker; /* marker array of each column in the panel */

                    /* Initialize pointers */
                    xsup       = Glu->xsup;
                    xsup_end   = Glu->xsup_end;
                    supno      = Glu->supno;
                    lsub       = Glu->lsub;
                    xlsub      = Glu->xlsub;
                    xlsub_end  = Glu->xlsub_end;
                    Astore     = (NCPformat *)A->Store;
                    a          = (doublecomplex*)Astore->nzval;
                    asub       = Astore->rowind;
                    xa_begin   = Astore->colbeg;
                    xa_end     = Astore->colend;
                    marker1    = marker + m;
                    repfnz_col = repfnz;
                    dense_col  = dense;
                    nextp      = 0;
                    *nseg      = 0;

                #if ( DEBUGlevel>=2 )
                if (jcol == BADPAN)    
                    printf("(%d) pzgstrf_panel_dfs[begin] jcol %d, w %d\n", pnum, jcol, w);
                #endif
    
                    /*
                     * For each column in the panel ...
                     */
                    for (jj = jcol; jj < jcol + w; ++jj, nextp += m) {
	                nextl_col = nextp;
	                col_marker = &spa_marker[nextp];

	                /*
	                 * For each nonz in A[*,jj] perform dfs ...
	                 */
	                for (k = xa_begin[jj]; k < xa_end[jj]; ++k) {
	                    krow = asub[k];
	                    dense_col[krow] = a[k];
	                    kmark = col_marker[krow];
	    
	                    /* if krow was visited before, go to the next nonzero */
	                    if ( kmark == jj ) continue;

	                    /*
	                     * For each unmarked nbr krow of jj ...
	                     */
	                    col_marker[krow] = jj;
	                    kperm = perm_r[krow];
	    
	                    if ( kperm == EMPTY ) {
		                /* krow is in L: place it in structure of L[*,jj].
		                 * NOTE: some entries in U may get here, because "perm_r"
		                 *       is not yet available from a preceeding busy column.
		                 */
		                panel_lsub[nextl_col++] = krow; /* krow is indexed into A */
	                    } else {
		                /* 
		                 * krow is in U (0 <= kperm < jcol): if its supernode
		                 * representative krep has been explored, update repfnz[*].
		                 */
		                if ( lbusy[kperm] == jcol ) { /* kperm is busy */
                #if ( DEBUGlevel>=3 )
                  if (jj == BADCOL)		    
                    printf("(%d) pzgstrf_panel_dfs(%d) skip busy krow %d, kperm %d\n",
	                   pnum, jj, krow, kperm);
                #endif		    
		                    continue;
		                }

		                /* Here, krep cannot possibly be "busy" */
		                krep = SUPER_REP( supno[kperm] );
		                myfnz = repfnz_col[krep];

                #ifdef CHK_DFS
                if (jj == BADCOL)		
                    printf("(%d) pzgstrf_panel_dfs[1] %d, krep %d, fsupc %d, Pr[krow %d] %d, myfnz %d\n",
	                   pnum, jj, krep, SUPER_FSUPC(supno[krep]), krow, kperm, myfnz);
                #endif
		                if ( myfnz != EMPTY ) {	/* Representative visited before */
		                    if ( myfnz > kperm ) repfnz_col[krep] = kperm;
		                    /* continue; */
		                } else {
		                    /* Otherwise, performs dfs starting from krep */
		                    parent[krep] = EMPTY;
		                    repfnz_col[krep] = kperm;
		                    if ( ispruned[krep] ) {
			                if ( SINGLETON( supno[krep] ) )
			                    xdfs = xlsub_end[krep];
			                else xdfs = xlsub[krep];
			                maxdfs = xprune[krep];
                #ifdef PROFILE
			                /*Gstat->procstat[pnum].pruned++;*/
                #endif		    
		                    } else {
			                fsupc = SUPER_FSUPC( supno[krep] );
			                xdfs = xlsub[fsupc] + krep-fsupc+1;
			                maxdfs = xlsub_end[fsupc];
                #ifdef PROFILE
			                /*Gstat->procstat[pnum].unpruned++;*/
                #endif		    
		                    }
                #ifdef CHK_DFS
                if (jj == BADCOL)		    
                {
                    register int i;
                    printf("(%d) pzgstrf_panel_dfs[2] %d, ispruned[%d] %d, xdfs %d, maxdfs %d\n",
	                   pnum, jj, krep, ispruned[krep], xdfs, maxdfs);
                    /*for (i = xdfs; i < maxdfs; i++) printf("(%d) lsub-%d", pnum, lsub[i]);*/
                    printf("\n");
                }
                #endif
		                    do {
			                while ( xdfs < maxdfs ) {
			                    /* for each unmarked kchild of krep ... */
			                    kchild = lsub[xdfs];
			                    xdfs++;
			                    chmark = col_marker[kchild];
			    
			                    if ( chmark != jj ) { /* Not reached yet */
				                col_marker[kchild] = jj;
				                chperm = perm_r[kchild];
				
				                if ( chperm == EMPTY ) {
				                    /* kchild is in L: place it in L[*,j]. */
				                    panel_lsub[nextl_col++] = kchild;
				                } else {
				                    /* kchild is in U (0 <= chperm < jcol): 
				                     * chrep = its supernode-rep. If its rep
				                     * has been explored, update its repfnz[*].
				                     */

				                    if ( lbusy[chperm] == jcol ) {
                #ifdef DEBUG
                if (jj == BADCOL)					
                    printf("(%d) pzgstrf_panel_dfs(%d) skip busy kchild %d, chperm %d\n",
	                   pnum, jj, kchild, chperm);
                #endif		    
	                                                     continue;
                                                    }
				    
				                    chrep = SUPER_REP( supno[chperm] );
				                    myfnz = repfnz_col[chrep];
                #ifdef DEBUG
                if (jj == BADCOL)				    
                    printf("(%d) pzgstrf_panel_dfs[3] %d, krep %d, Pr[kchild %d] %d, chrep %d, fsupc %d, myfnz %d\n",
	                   pnum, jj, krep, kchild, chperm, chrep,
	                   SUPER_FSUPC(supno[chrep]), myfnz);
                #endif
				                    if ( myfnz != EMPTY ) {/* Visited before */
					                if ( myfnz > chperm )
					                    repfnz_col[chrep] = chperm;
				                    } else {
					                /* Cont. dfs at snode-rep of kchild */
					                xplore[krep] = xdfs;
					                xplore[m + krep] = maxdfs;
					                parent[chrep] = krep;
					                krep = chrep; /* Go deeper down G(L) */
					                repfnz_col[krep] = chperm;
					                if ( ispruned[krep] ) {
					                    if ( SINGLETON( supno[krep] ) )
						                xdfs = xlsub_end[krep];
					                    else xdfs = xlsub[krep];
					                    maxdfs = xprune[krep];
                #ifdef PROFILE
					                    /*procstat[pnum].pruned++;*/
                #endif		    
					                } else {
					                    fsupc = SUPER_FSUPC(supno[krep]);
					                    xdfs = xlsub[fsupc] + krep-fsupc+1;
					                    maxdfs = xlsub_end[fsupc];
                #ifdef PROFILE
					                    /*procstat[pnum].unpruned++;*/
                #endif		    
					                }
                #ifdef CHK_DFS
                if (jj == BADCOL)
                    printf("(%d) pzgstrf_panel_dfs[4] %d, ispruned[%d] %d, xdfs %d, maxdfs %d\n",
	                   pnum, jj, krep, ispruned[krep], xdfs, maxdfs);
                #endif
					
				                    } /* else */
				                } /* else */
			      
			                    } /* if... */
			    
			                } /* while xdfs < maxdfs */
			
			                /* krow has no more unexplored nbrs:
			                 *    Place snode-rep krep in postorder DFS, if this 
			                 *    segment is seen for the first time. (Note that
			                 *    "repfnz[krep]" may change later.)
			                 *    Backtrack dfs to its parent.
			                 */
			                if ( marker1[krep] != jcol ) {
			                    segrep[*nseg] = krep;
			                    ++(*nseg);
			                    marker1[krep] = jcol;
                #ifdef CHK_DFS
                if (jj == BADCOL)			    
                    printf("(%d) pzgstrf_panel_dfs(%d) repfnz[%d] %d added to top.list by jj %d\n",
	                   pnum, jj, krep, repfnz_col[krep], jj);
                #endif			    
			                }
			
			                kpar = parent[krep]; /* Pop stack, mimic recursion */
			                if ( kpar == EMPTY ) break; /* dfs done */
			                krep = kpar;
			                xdfs = xplore[krep];
			                maxdfs = xplore[m + krep];
			
                #ifdef CHK_DFS
                if (jj == BADCOL)			
                {
                    register int i;
                    printf("(%d) pzgstrf_panel_dfs[5] pop stack: %d, krep %d, xdfs %d, maxdfs %d\n",
	                   pnum, jj, krep, xdfs, maxdfs);
                    /* for (i = xdfs; i < maxdfs; i++) printf("(%d) lsub-%d", pnum, lsub[i]);*/
                    printf("\n");
                }
                #endif

		                    } while ( kpar != EMPTY ); /* until empty stack */
		    
		                } /* else: myfnz == EMPTY */
		
	                    } /* else: kperm != EMPTY */
	    
	                } /* for each nonzero in A[*,jj] */

                #if ( DEBUGlevel>=3 )
                if (jj == BADCOL) {
                #define REPCOL 0    
                    krep = REPCOL;
                    printf("(%d) pzgstrf_panel_dfs(end) w_lsub_end[jj=%d] %d, repfnz_col[%d] %d\n",
	                   pnum, jj, nextl_col - nextp, krep, repfnz_col[krep]);
                    PrintInt10("lsub", nextl_col - nextp, &panel_lsub[nextp]);
                }
                #endif
	
	                w_lsub_end[jj-jcol] = nextl_col - nextp;
	                repfnz_col += m;
                        dense_col += m;
	
                    } /* for jj ... */

                }

            };
        };
    };
};
